import numpy as np
import pandas as pd
import plotly.graph_objects as go
from statsmodels.tsa.statespace.sarimax import SARIMAX
import warnings
# 警告を無視（収束に関する警告などが出る場合があるため）
warnings.filterwarnings("ignore")SARIMAを使った医療費予測のデモ
Statistics
Time Series
Demo
SARIMA（Seasonal AutoRegressive Integrated Moving Average）は、季節性のある時系列データに適した予測モデルです。医療費の階差を取ると、定常分布に近づき、医療費は冬はインフルエンザが流行したり夏は熱中症が増えたりと季節性があるため、SARIMAモデルを使用して医療費の予測を行います。
1. データの準備
# 1. データの作成
np.random.seed(42)
# 期間: 過去10年分 (120ヶ月)
dates = pd.date_range(start='2015-01-01', periods=120, freq='MS')
# トレンド: 徐々に医療費が上がっていく (線形トレンド)
trend = np.linspace(100, 150, 120) # 100単位から1000単位へ上昇 <- ここを変えた。
# 季節性: 12ヶ月周期 (例: 冬に高く、夏に低いなど)
seasonality = 10 * np.sin(np.linspace(0, 20 * np.pi, 120))
# ノイズ: ランダムな変動
noise = np.random.normal(scale=3, size=120)
# 合成して「月次医療費データ」とする
medical_costs = trend + seasonality + noise
# DataFrame化
df = pd.DataFrame({'Date': dates, 'Cost (100 million)': medical_costs})
df.set_index('Date', inplace=True)
df["Cost (100 million)"] = df["Cost (100 million)"].round(2)
print("データの先頭5行:")
display(df.head())データの先頭5行:| Cost (100 million) | |
|---|---|
| Date | |
| 2015-01-01 | 101.49 |
| 2015-02-01 | 187.82 |
| 2015-03-01 | 277.03 |
| 2015-04-01 | 364.15 |
| 2015-05-01 | 440.64 |
実績データの可視化 (Plotly)
# 2. 実績データの可視化
fig = go.Figure()
fig.add_trace(go.Scatter(
x=df.index,
y=df['Cost (100 million)'],
mode='lines+markers',
name='実績値',
line=dict(color='blue')
))
fig.update_layout(
title='医療費ダミーデータ（過去10年）',
xaxis_title='年月',
yaxis_title='医療費（1億円）',
template='plotly_white'
)
fig.show(renderer="notebook")SARIMAモデルの構築と学習
\(\phi_p(B) \Phi_P(B^s) (1 - B)^d (1 - B^s)^D y_t = \theta_q(B) \Theta_Q(B^s) \varepsilon_t\)
SARIMAモデル（季節性自己回帰和分移動平均モデル）を構築します。
※本来はAICなどでパラメータ探索を行いますが、今回はダミーデータの構造に合わせて一般的なパラメータ \((p, d, q) \times (P, D, Q, s)\) を設定します。
# 3. SARIMAモデルの構築
# order=(p, d, q), seasonal_order=(P, D, Q, s)
# 周期 s=12 (月次データのため)
sarima_model = SARIMAX(
df['Cost (100 million)'],
order=(1, 1, 1),
seasonal_order=(1, 1, 1, 12),
enforce_stationarity=False,
enforce_invertibility=False
)
# モデルの学習
sarima_result = sarima_model.fit()
print(sarima_result.summary()) SARIMAX Results
==========================================================================================
Dep. Variable: Cost (100 million) No. Observations: 120
Model: SARIMAX(1, 1, 1)x(1, 1, 1, 12) Log Likelihood -244.124
Date: Sun, 14 Dec 2025 AIC 498.249
Time: 15:21:45 BIC 510.912
Sample: 01-01-2015 HQIC 503.362
- 12-01-2024
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.0269 0.121 0.223 0.824 -0.210 0.263
ma.L1 -1.0000 238.268 -0.004 0.997 -467.996 465.996
ar.S.L12 -0.2240 0.144 -1.554 0.120 -0.507 0.058
ma.S.L12 -0.5145 0.166 -3.101 0.002 -0.840 -0.189
sigma2 10.2597 2444.753 0.004 0.997 -4781.367 4801.887
===================================================================================
Ljung-Box (L1) (Q): 0.01 Jarque-Bera (JB): 0.43
Prob(Q): 0.91 Prob(JB): 0.81
Heteroskedasticity (H): 0.95 Skew: 0.08
Prob(H) (two-sided): 0.88 Kurtosis: 3.29
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).向こう5年間の予測と可視化
# 4. 向こう5年間 (60ヶ月) の予測
forecast_steps = 60
pred_uc = sarima_result.get_forecast(steps=forecast_steps)
# 予測値（期待値）
pred_mean = pred_uc.predicted_mean
# 95%信頼区間
pred_ci = pred_uc.conf_int()
# --- 可視化 ---
fig_forecast = go.Figure()
# 1. 実績値のプロット
fig_forecast.add_trace(go.Scatter(
x=df.index,
y=df['Cost (100 million)'],
mode='lines',
name='実績値',
line=dict(color='blue')
))
# 2. 予測値のプロット
fig_forecast.add_trace(go.Scatter(
x=pred_mean.index,
y=pred_mean,
mode='lines',
name='予測値 (今後5年)',
line=dict(color='red', dash='dash')
))
# 3. 信頼区間のプロット (上限と下限の間を塗りつぶす)
fig_forecast.add_trace(go.Scatter(
x=pred_ci.index,
y=pred_ci.iloc[:, 0], # 下限
mode='lines',
line=dict(width=0),
showlegend=False,
name='Lower Bound'
))
fig_forecast.add_trace(go.Scatter(
x=pred_ci.index,
y=pred_ci.iloc[:, 1], # 上限
mode='lines',
line=dict(width=0),
fill='tonexty', # ひとつ前のトレース（下限）との間を埋める
fillcolor='rgba(255, 0, 0, 0.2)', # 赤色の半透明
name='95%信頼区間'
))
fig_forecast.update_layout(
title='SARIMAモデルによる医療費予測（向こう5年）',
xaxis_title='年月',
yaxis_title='医療費',
template='plotly_white',
hovermode="x unified"
)
fig_forecast.show()Unable to display output for mime type(s): application/vnd.plotly.v1+json医療費上昇金額の算出
# 5. 上昇金額の算出
# 直近1年間の実績合計
last_year_actual = df['Cost (100 million)'].iloc[-12:].sum()
# 5年後（予測期間の最後の1年間）の予測合計
last_year_forecast = pred_mean.iloc[-12:].sum()
# 上昇額と上昇率
increase_amount = last_year_forecast - last_year_actual
increase_rate = (increase_amount / last_year_actual) * 100
print(f"--- 医療費上昇シミュレーション ---")
print(f"直近1年間の医療費合計: {last_year_actual:,.2f}")
print(f"5年後の年間医療費合計: {last_year_forecast:,.2f}")
print(f"--------------------------------")
print(f"年間増加額: +{increase_amount:,.2f}")
print(f"上昇率: {increase_rate:.2f}%")--- 医療費上昇シミュレーション ---
直近1年間の医療費合計: 1,777.19
5年後の年間医療費合計: 2,088.13
--------------------------------
年間増加額: +310.94
上昇率: 17.50%